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ABSTRACT 

The torques exerted by a locally isothermal disk on an embedded planet lead to rapid inward migration. Recent 
work has shown that modeling the thermodynamics without the assumption of local isothermality reveals regions 
where the net torque on an embedded planet is positive, leading to outward migration of the planet. When a region 
with negative torque lies directly exterior to this, planets in the inner region migrate outwards and planets in the 
outer region migrate inwards, converging where the torque is zero. We incorporate the torques from an evolving 
non-isothermal disk into an N-body simulation to examine the behavior of planets or planetary embryos interacting 
in the convergence zone. We find that mutual interactions do not eject objects from the convergence zone. Small 
numbers of objects in a laminar disk settle into near resonant orbits that remain stable over the 10 Myr periods that 
we examine. However, either or both increasing the number of planets or including a correlated, stochastic force to 
represent turbulence drives orbit crossings and mergers in the convergence zone. These processes can build gas giant 
cores with masses of order ten Earth masses from sub-Earth mass embryos in 2-3 Myr. 



1. INTRODUCTION 

The details of how giant planets form in protoplan- 
etary disks remain poorly understood. The core ac- 
cretion model developed by Pollack et al. (1996) re- 
lies on the rapid growth of a gas giant through run- 
away accretion of gas onto a solid planetary core that 
grows massive enough to gravitationally capture the 
gas. Such cores form through binary collisions of plan- 
etary embryos — the bodies formed during oligarchic 
growth, with masses ranging from 0.1 to a few Earth 
masses. Since the accretion must necessarily occur in 
a gas-rich disk, the planetary core must be formed be- 
fore the gas disk dissipates. The mass required to begin 
rapidly accreting gas is determined by the aspect ratio 
of the disk, the mass of the star, and the luminosity of 
the core that is provided by planetesimal accretion and 
can prevent collapse of the accreting gas. For a solar 
mass star, a thin disk with an aspect ratio of 0.05, and 
a surface density of planetesimals of 10 g cm -2 at 5 AU, 
the core must be around ten Earth masses to begin run- 
away accretion. The lifetime of the gas disk is on the 
order of 10 Myr, and the core must form before the disk 
dissipates. 

Planetary embryos are massive enough to not be af- 
fected by the gas drag forces that dominate the motion 
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of smaller bodies moving through the sub-Keplerian gas 
disk. The mass of the planetary embryos, however, in- 
duces spiral Lindblad resonances in the gas disk. The 
asymmetry in the size of the inner and outer Lindblad 
resonances generates a net negative torque from the 
disk that, in the absence of other effects, robs the planet 
of angular momentum, leading to inward Type I migra- 
tion and ultimately infall into the central star. A linear 
treatment of the torque from isothermal gas in the co- 
rotating region of the planetary embryo shows that it 
provides a positive torque arising from the presence of 
a vortensity gradient that nearly offsets the net Lindblad 
torque, but a small net negative torque remains. 

The behavior of multiple planetary embryos embed- 
ded in a gas disk and undergoing inward migration 
due to Lindblad and linear co-rotational torques was ex- 
plored by Cresswell & Nelson (2008). They found that 
as bodies with different masses migrate inwards, differ- 
ential migration leads to mergers, forming bodies with 
masses in the planetary core range. These simulations 
use large embryos-with masses of at least 2 M© and up 
to planetary core mass-and simulate migration in a lam- 
inar, isothermal disk. 

In addition to the perturbations to the gas disk caused 
by the presence of the planetary embryos, a partially 
ionized disk with a subthermal magnetic field can de- 
velop turbulence through the magnetorotational insta- 
bility (MRI, Balbus & Hawley 1991). The resulting ran- 
dom overdensities in the gas can lead the planetary 
embryo to undergo a random walk, possibly avoid- 
ing infall for some of the bodies (Nelson 2005). While 
this is an effective way to prevent catastrophic infall 
for planetesimal-sized objects with radii of < 1000 km 
(Yang et al. 2009, 2011), it does not provide a clear mech- 
anism for building large planetary cores. 

Planetary embryos undergoing inward Type I migra- 
tion can be trapped in pressure ridges that may form in 
the presence of a dead-zone boundary or within an anti- 
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cyclonic vortex. The positive gradient in surface den- 
sity in these regions leads to formation of a convergence 
zone, where migration halts. The dynamics of multiple 
planetary embryos undergoing Type I migration in the 
proximity of a pressure ridge convergence zone have 
been investigated by Morbidelli et al. (2008), who found 
that planets collide until only a few remain. Sandor 
et al. (2011) showed that such pressure traps can lead 
to growth from planetary embryos to planetary cores 
through continual generation of bodies that migrate into 
the pressure ridge convergence zone. 

One possible means of avoiding catastrophic infall 
that could lead to rapid formation of planetary cores 
was found by Paardekooper & Mellema (2006), who ex- 
plored the interactions of planets with a gas disk, in a 
model in which they did not make the locally isothermal 
approximation, but rather directly modeled the effects 
of radiative transfer through optically thick gas with 
an adiabatic equation of state. They found that the co- 
rotation torque arising from the gas undergoing horse- 
shoe turns near the planet showed a strong positive con- 
tribution when the isothermal approximation was re- 
laxed. This was later shown to be a result of the pres- 
ence of an entropy gradient across the corotation region 
(Paardekooper & Papaloizou 2008; Baruteau & Masset 
2008). A large radial entropy gradient increases the pos- 
itive co-rotation torque to the point that it cancels out 
the negative net Lindblad torque and can even reverse 
the direction of migration. This understanding was de- 
veloped further in Paardekooper & Papaloizou (2009). 
A prescription for the nonlinear horseshoe torque was 
given by Paardekooper et al. (2010), allowing the to- 
tal torque on a planet undergoing Type I migration to 
be calculated, given the local temperature and surface 
density gradients. The migration of a single planet em- 
bedded in an evolving disk was then studied by Lyra et 
al. (2010, hereafter LPM10), who found that planets mi- 
grate towards convergence zones, where the net torque 
is zero. Planets located exterior to these orbits experi- 
ence a negative torque, while planets located interior 
experience a positive torque. Upon reaching the zero- 
torque region, the planets migrate inwards with the disk 
until the disk dissipates, thus avoiding collision with the 
star. 

In this paper, we extend the work of LPM10 by fol- 
lowing the interactions of multiple planets converging 
towards any zero-torque radius. We use the disk evo- 
lution model from LPM10, for a non-isothermal one- 
dimensional disk with photoevaporation and viscous 
diffusion. The radial profiles of temperature and surface 
density from the disk evolution model are used to de- 
termine the torque on planetary embryos in an N-body 
simulation, following the formula for non-isothermal, 
unsaturated, horseshoe torque of Paardekooper et al. 
(2010). An optimized version of the Bulirsch-Stoer N- 
body code used in Sandor et al. (2011) is employed to 
integrate the orbits of the planets under the influence 
of the torque from the gas disk, as well as eccentricity 
and inclination damping following Cresswell & Nelson 
(2008). The code was modified to include the effect of 
turbulence on the planetary motion, using the prescrip- 
tion described by Laughlin et al. (1994). 

We start in Sect. 2 by describing the disk model, the 
prescriptions for the torque, eccentricity and inclination 
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FIG. 1. — The dimensionless migration torque and how it evolves 
through time. The region of interest in these simulations is the tran- 
sition from negative torque - inward migration - to positive torque 
that is initially located around 22 AU. As the disk evolves, this con- 
vergence zone moves inwards. The torque shown in this figure is an 
interpolation between the adiabatic and isothermal torques, based on 
the opacity at a given radius. To find the torque acting on a body at 
a given location, this dimensionless torque must be multiplied by To, 
from Eq. (12). 

damping, the turbulence model, and how they are in- 
cluded in the N-body code. In Sect. 3, we describe the 
initial conditions of the planets in the simulations as 
well as the parameters describing the different gas disks 
that were used. Section 4 presents the results of the sim- 
ulations and Sect. 5 includes a discussion of the results. 
A study tackling a similar problem was made public af- 
ter initial submission of our paper by Hellary & Nelson 
(2012, hereafter HN12). We discuss the differences be- 
tween our work and theirs in Sect. 5.1. A summary and 
conclusions are presented in Sect. 6. 

2. METHOD 

Our simulation can be viewed as a modified N-body 
simulation, which incorporates the migration torques 
acting on planetary embryos, a stochastic force repre- 
senting turbulence in the gas disk, and eccentricity and 
inclination damping from the gas disk, as well as the 
disk evolution due to viscosity and photoevaporation. 
The migration torques depend on temperature and sur- 
face density profiles from one-dimensional simulations 
of the evolution of the gas disk, and the force from tur- 
bulence is modeled using a time-dependent potential. 
Each component of the simulation is described in more 
detail in this section. 

2.1. Disk Evolution 

The gas disk evolves via viscous diffusion and pho- 
toevaporation. The evolution of the gas surface density 
due to viscosity is given by (Lynden-Bell & Pringle 1974) 



~dt 



3d_ 
r dr 



A/2 



dZvr 1/2 
dr 



(1) 



where r is the distance from the central star, and v is 
the kinematic viscosity. Following Shakura & Sunyaev 
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TABLE 1 

Run Parameters 



Run 


M„ 


f mass 


*^mass 


Combine 


N p 


N m HR 


Mtot 


Swarm 


Tform 




(1) 


PJ 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(ii) 


HI (Fiducial) 


0.75 


- 


- 


- 


23 


3 


17.25 


No 


1.29 


10.5 


H2 


0.5 


— 


— 


— 


27 


3 


13.5 


No 


2.23 


8.0 


H3 


1.00 


— 


— 


— 


21 


3 


21.0 


No 


0.43 


13.0 


H4 


1.25 


— 


— 


— 


20 


3 


25.0 


No 


0.21 


10.0 


H5 


1.5 


— 


— 


— 


19 


3 


29.5 


No 


0.15 


10.5 


Gl 


— 


1.0 


0.3 


— 


21 


3 


20.9 


No 


0.27 


11.1 


bZ 


- 


1.0 


0.3 




11 


6 


11.7 


No 


0.34 


7.1 


G3 


— 


0.75 


0.3 


— 


24 


3 


18.7 


No 


2.38 


10.2 


G4 


— 


0.75 


0.3 


— 


13 


6 


9.8 


No 


1.49 


9.3 


G5 




0.5 


0.2 




27 


3 


13.9 


No 


3.18 


9.3 


G6 




0.5 


0.2 




14 


6 


6.2 


No 


1.38 


4.4 


G7 








G2+G6 


2b 


6 


17.9 


No 


1.36 


12.0 


G8 








G2+G4 


24 


6 


21.5 


No 


0.95 


10.5 


G9 








G4+G6 


27 


6 


16.0 


No 


2.64 


9.4 


Swl 


0.25 








25 


4 


6.25 


Yes 


4.60 


2.5 


Sw2 


0.75 








12 


6 


9 


Yes 


4.41 


8.5 


Sw3 


1.0 








11 


6 


11 


Yes 


0.68 


8.5 


Sw4 


1.0 








7 


10 


7 


Yes 


1.42 


9.0 


SD1 


0.75 








23 




17.25 


No 


0.98 


10.5 


SD2 


1.0 








17 




17.0 


No 


0.267 


7.0 & 5.0 


SD3 




0.75 


0.3 




21 




17.26 


No 


1.03 


10.2 



NOTE. — Col. (1): Name of run. Col. (2): Initial mass of bodies, in . Col. (3): Mean of Gaussian distribution for initial masses, in . Col. (4): Standard 
deviation of Gaussian distribution for initial masses, in M^. Col. (5): Combined initial sample. Col. (6): Number of initial planets. Col. (7): initial separation 
between planets, in mutual Hill radii. Col. (8): Total initial mass. Col. (9): Presence of background swarm of small (0.25 Mj) planets. Col. (10): Formation time 
(Myr). Col. (11): Mass of the core formed, in M@ . 



(1973) we use 



v — aci/fl 



(2) 



where a is a dimensionless parameter. The Keplerian 
angular velocity is 



n = (GM*/r 



3x1/2 



(3) 



where is the mass of the central star, and c s is the lo- 
cal sound speed. Our simulation uses a model based on 
the results of Hollenbach et al. (1994) to specify the rate 
of surface density loss due to photoevaporation. This 
occurs only beyond a radius r g where ionized disk gas 
has a sound speed equal to the escape velocity. In this 
regime 



M w rl 5 



2n(r ext -r g )r 2 - 5 ' 



(4) 



while S w — within. The rate of photoevaporation is 
set by the wind mass loss rate M w , the escape radius 
r g , and the external boundary of our simulation r ext . As 
the gas disk evolves, we record temperature and surface 
density profiles every 10 3 yr. 

For temperature evolution, we use a model without 
shock heating (Nakamoto & Nakagawa 1994) that bal- 
ances viscous heating, background radiation, and radia- 
tive cooling to derive an emerging flux 



T eff 



9 9 

4 



+ 2*1*, 



(5) 



where T and T fc are the midplane and background tem- 
peratures, and a the Stefan-Boltzmann constant. We 
take the effective optical depth at the midplane (Hubeny 



1990; Kley & Crida 2008) 



3t y/3 1 

T ^ = y + X + 4T 



(6) 



The optical depth is r = kS/2, and the opacities k 
are taken from Bell et al. (1997). We assume that al- 
though dust growth and planet formation lock away re- 
fractory material, fragmentation efficiently replenishes 
small grains, keeping the disks opaque during their evo- 
lution (Birnstiel et al. 2009). 

The gas simulations are conducted using the same nu- 
merical algorithm as LPM10, which evolves the gas disk 
in one dimension over the lifetime of the disk. As a nu- 
merical solver, we used the PENCIL CODE 7 , which inte- 
grates the equations of evolution with sixth order spatial 
derivatives, and a third order Runge-Kutta time integra- 
tor. Boundary conditions are taken as outflow. Because 
the optical depth depends on temperature, we solve Eq. 
(5) with a Newton-Raphson root-finding algorithm (us- 
ing 0.01 K precision). To save on computational time, 
temperatures are pre-computed as a function of £ and 
Q, stored in the memory as look-up tables, and retrieved 
in run time via bilinear interpolation. 

The reduction to one dimension assumes azimuthal 
symmetry and that the disk can be vertically integrated 
without significantly changing the characteristics of the 
disk. The azimuthal asymmetry will be reintroduced 
into the simulation via our migration and turbulence 
prescriptions, and the scale-height parameterizes the ef- 
fects of a finite disk height on the magnitude of the tur- 
bulence as well as setting the timescales for eccentricity 
and inclination damping. 

7 The code, including improvements done for the present work, is 
publicly available under a GNU open source license and can be down- 
loaded at http:/ / www.nordita.org/ software/ pencil-code 
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The profiles from the gas simulations are used to de- 
termine the local temperature and gas surface density in 
the N-body code. The migration torques and turbulent 
forces acting on the planets are proportional to the local 
surface density E and temperature T. The local values 
of E and T at any given time and position are found 
by linearly interpolating from the two closest times, and 
then again between the two nearest radial grid points. 

Several quantities are indirectly dependent on T by 
way of the local sound speed, 



(7) 



where 7 is the adiabatic coefficient of the gas, set to 7/5 
for diatomic molecules, k is Boltzmann's constant, ]i is 
the mean molecular weight of the gas, set to 2.4 for a 5:1 
H2-He mixture, and OTh is the mass of a hydrogen atom. 
The sound speed is also used to determine the aspect 
ratio of the disk 

r 




vk 



(8) 



where H is the scale height of the disk, and the Keple- 
rian velocity v% = (GM^/r) 1 ^ 2 . Unless otherwise noted, 
E, c s , H and h are time and radius dependent, with val- 
ues derived from the gas evolution model. 

The local radial gradients of temperature and surface 
density needed for the torque calculation are calculated 
for the entire disk during its evolution using a sixth- 
order, finite-difference method, as is done in the PEN- 
CIL CODE. The gradients are recorded with the density 
and temperature profiles. During the N-body simula- 
tion, the same linear interpolation method is used for 
the gradients as for the original quantities. 

We ignore gap opening by massive planets, since that 
does not occur for planets with masses less than 10- 
20 M© (Crida et al. 2006), and the goal of this simula- 
tion is to examine whether planetary embryos can coa- 
lesce into planets with masses approaching this regime. 
Therefore, we can simplify the simulations by evolving 
the gas disk for its entire lifetime before beginning the 
N-body simulations. 

2.2. Prescription for migration 

Planets in a gas disk induce spiral Lindblad reso- 
nances in the disk, as well as driving the gas in the co- 
orbital region into horseshoe orbits. The net torque on 
the planet from the Lindblad resonances can be calcu- 
lated analytically, as in Ward (1997). In isothermal, non- 
magnetized disks, the full net torque (Lindblad and co- 
rotational) can be calculated semi-analytically (Tanaka 
& Ward 2004). The Lindblad resonances impart a neg- 
ative torque to the planet, leading to the well known 
effect of rapid inward migration. Hydrodynamical sim- 
ulations of a planet with mass on the order of M© em- 
bedded in an isothermal disk have confirmed the pre- 
dictions from linear theory (Bate et al. 2003; D'Angelo et 
al. 2002). 

Paardekooper & Mellema (2006) investigated the mi- 
gration of a low-mass planet in an adiabatic, hydrody- 
namic simulation including radiative transfer for cool- 
ing radiation, and showed that outward migration can 
occur in such a non-isothermal disk due to positive 



torques exerted by the co-rotation region under some 
circumstances. Later investigations showed that this de- 
parture from isothermal Type I migration is caused by 
the radial entropy gradient of the disk (Paardekooper 
& Papaloizou 2008; Baruteau & Masset 2008). A pre- 
scription incorporating the effects of non-isothermal co- 
rotation torques was developed by Paardekooper et al. 
(2010). 

The prescription allows us to calculate the torque act- 
ing on a low-mass planet embedded in a gas disk as a 
function of the radial gradient of the disk's temperature 
and density in the region near the planet. Defining the 
gradients as 

dlnE „ dlnT 

1= ain7 ; ^ = ^ ; £ = P-(7-i)^ (9) 



<91nr ' 



Paardekooper et al. (2010) arrived at the following equa- 
tion for the total unsaturated torque on the planet 



7r ad /ro = -0.85 -n- 1.75/3 + 7.%/ 7 . 



(10) 



If the gas passing through the horseshoe turn cools 
rapidly enough, the co-rotational torque is best approx- 
imated by the isothermal model 

r iS o/r =-0.85-/7-0.9/5. (11) 

The scale of both torques is set by 

r =(i)W, 



(12) 



where q is the ratio of the planet mass to the central stel- 
lar mass, is the distance from the planet to the central 
body. This prescription is valid for low mass planets in 
unsaturated regions. 

If the gas can radiate away the surplus heat and re- 
store the temperature quickly, then the gas exerting the 
horseshoe drag should be treated isothermally. The ef- 
fective torque is interpolated between the isothermal 
and adiabatic torque models. The ratio of the radiative 
cooling time, t md , to the dynamical time, t dyn , defines 
the interpolation between the two regimes. We estimate 



trad ~E/E, 



(13) 



where £ is the rate of cooling by radiation, or equiva- 
lently the divergence of the radiation flux,F, and E is 
the internal energy density of the gas, E = cypT. For a 
cylinder of height 2H centered on the disk, p = E/2H, 
and the total energy within the cylinder is 



EdV - 



CyET 
2H 



dV = nr CyET, 



(14) 



where the last step is true for a cylinder with a small 
radius with constant temperature and surface density 
throughout. Integrating the rate of cooling over the 
same volume gives 



J EdV = J V • FdV = ^FdS : 



-2nr 2 o-T^. 



(15) 



This is true for small radii if the temperature outside of 
the cylinder is nearly the same as the temperature in- 
side, since then the flux is approximately zero every- 
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where except the top and bottom of the cylinder. Sub- 
stituting T 4 ff = T 4 / T e ff, the radiation timescale is 



£ 



CyST eff 



(16) 



Defining to be the ratio of the radiative timescale to 
the dynamic timescale 2k/ Q, we have 



= 



CyZ']7T e ff 

AnaT 3 



(17) 



The effective torque acting on a planet is then given 
by interpolating between the adiabatic and isothermal 
torques 

Tad© 2 + Tiso 



typel 



(0 + 1)' 



(18) 



This interpolation, though not derived from physical 
principles, provides a smooth mathematical transition 
between the isothermal and adiabatic limits. We stress, 
though, that the precise nature of the transition has lit- 
tle influence on the results since, as shown in LPM10, 
most of the planet evolution in the disk through space 
and time occurs well within the adiabatic regime. The 
resulting torques are shown in Fig. 1. 

2.3. Prescription for eccentricity and inclination damping 

In addition to the torque that is responsible for Type I 
migration, the gas disk also imparts a force on a planet 
that acts to dampen the planet's eccentricity e and in- 
clination i (Tanaka & Ward 2004). This dampening can 
be modeled as an exponential decay. Tanaka and Ward 
(2004) provide the timescale 



Mih 



2 1,4 



amp 



nipSa^Sl 



(19) 



where a-p is the semi-major axis of the planet, and the 
disk aspect ratio h is given by equation (8). Cresswell 
& Nelson (2008) showed that the exponential decay is 
only applicable for small inclinations and eccentricities, 
after which point the decay is best modeled as a power 
law. Following their prescriptions, we use the following 
timescales for the dampening of eccentricity and incli- 
nation, respectively 



te = ^ ( l _ .14£ 2 + 0.06e 3 + 0.18« 2 



0.780 



(20) 



1 - 0.30 i 2 + 0.24 1 3 + 0.14 is 2 



(21) 



^damp 

! ' ~ 0.544 

where e=e/h and t=i/h. 

These equations give a timescale for the damping that 
agrees well with hydrodynamic simulations over the 
range of eccentricities and inclinations observed in our 
simulations 8 . The timescales are calculated at each step 

8 These timescales were derived for isothermal disks. However, 
Bitsch & Kley (2010, 2011) show that the behavior in radiative disks 
is qualitatively similar. The quantitative difference is that the damp- 
ing timescale is slightly longer in radiative disks, because of the higher 
sound speed in an adiabatic gas. This is expected, since the damping 
is provided by backreaction of waves generated in the disk by an ec- 
centric and /or inclined planet. 




FIG. 2. — The history of our fiducial model (HI), initially consisting 
of 23 planetary embryos of 0.75M^ each. The planets migrate towards 
the convergence zone. Chaotic N-body interactions then start, and 
continue until a massive core is formed. The planetary core eventu- 
ally reaches a mass of 13M® and continues to occupy the zero-torque 
region, with smaller bodies interior and exterior to it being forced into 
mean-morion resonances. 
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FIG. 3. — The fiducial model, (HI), creates a large core that eventu- 
ally reaches a mass of 13Me- The remaining bodies are trapped in 
resonances with the planetary core. Following the transition at 1.57 
Myr, the resonances remain stable throughout the rest of the simula- 
tion, with two 3Mj, bodies trapped in the 5:6 resonance, a 6M^ body 
in a 7:6 resonance, and a 3M^ embryo trapped as a Trojan of the plan- 
etary core. Similar results are seen in other runs. 

of the N-body integration, and the force from the damp- 
ing effect on a planet moving with velocity v is calcu- 
lated as 



(v-r)r 

* damp,r = — 2 ^ MpV 



damp,z , TltpZ, 



(22) 
(23) 



where r and z are the unit vectors in the r and z direc- 
tions, respectively. 

2.4. Turbulence 
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FIG. 4. — The time required to form a core of at least 10 Me is shown 
to be a function of the initial masses of the planetary embryos. These 
runs use a uniform mass for all bodies, whose value sets the x coordi- 
nate. All initial conditions except those for the 1.5M® embryos have 
one simulation of the 5 runs fail to form a 10 M® core. These models 
are not included in the figure. 

Small enough planets embedded in a turbulent disk 
have been shown to undergo random walks, rather than 
migrate inward in classical Type I migration (Nelson & 
Papaloizou 2004; Nelson 2005). The random walks can 
be modeled as a first order Markov process (Rein & Pa- 
paloizou 2009) where the root mean squared torque and 
the autocorrelation time can be calculated up to a con- 
stant by dimensional analysis. However, the correlation 
function of such a first order Markov process over the 
plane of the disk will be zero for all non-zero distances — 
so two planets in close proximity may feel completely 
different forces from this prescription for turbulence. 

We use a model developed by Laughlin et al. (1994) 
and further modified by Ogihara et al. (2007) that treats 
the turbulent force as the gradient of a potential mul- 
tiplied by a scaling factor. In this model, the turbulent 
force is continuous over the disk and correlated torques 
act on planets near to each other. Due to the chaotic na- 
ture of N-body problems, turbulence can lead to colli- 
sions that would not have occurred otherwise, but dur- 
ing the close encounter of two bodies, the difference 
between the turbulent forces acting on each body goes 
to zero, and the relative velocity at collision is close to 
the velocity for two bodies on the same orbit in the ab- 
sence of turbulence. Additionally, our results show that 
neighboring orbits near the convergence zone can be 
stable when separated by as little as 0.4 AU at a distance 
of 11 AU from the central star, so the same turbulent per- 
turbations could affect both orbits. 

The turbulent potential, <P, is described as the sum of a 
large number of independent, scaled oscillation modes 



ipr 2 f2 2 A C/n 



(24) 



where A c m is one oscillation mode, defined as 



A c , m = ^g-(^) /(7 cos (m6 -<p c - Q c t) sin (^-j (25) 
Each mode is defined by the azimuthal wavenumber 



f 



m, as well as r c and (p c , which specify the initial cen- 
ter of the mode. The magnitude of the perturbation is 
set by £, whose value is set by a normal distribution 
with standard deviation of 1 and a mean value of 0. The 
planet's location is given by the radial coordinates r and 
6, and the displacement above the disk, z, is assumed to 
be small enough to have a negligible effect. The mode 
comes into existence at time tg and evolves as a func- 
tion of t = to + t. The lifetime of the perturbation is 
At — 2nr c / (mc s ), which is the sound-crossing time for a 
mode, where c s is the speed of sound. In Eq. (25), fl c is 
the Keplerian angular velocity at r c . 

The value for m is chosen from a log-random distribu- 
tion ranging from 1 to 64. The value of r c is selected from 
a uniform random distribution ranging from r- mt to r ext , 
and (p c is selected from a uniform distribution from to 
In. The radial profile of the perturbation is a Gaussian 
whose size is set by a = nr c f (Am). At the beginning of 
a simulation, n independent modes are generated, and 
as one mode expires another is generated so that at any 
time there are n perturbations. This parameter n is set 
to 50 in Laughlin et al. (1994) as well as in Ogihara et 
al. (2007), but since we are investigating a disk that is 
an order of magnitude larger than that in Laughlin et al. 
(1994), we have included 100 modes. 

In Eq. (24), the magnitude of the potential used to 
model turbulence is determined by the dimensionless 
parameter ip. This sets the scale for the magnitude of the 
turbulent force in comparison to the laminar forces due 
to Lindblad resonances and horseshoe torques. Ogi- 
hara et al. (2007) suggest that ip can be anywhere from 
10~ 3 to 10" 1 , depending on the strength of the instability 
driving the turbulence. The relation between ip and the 
Shakura & Sunyaev (1973) a parameter was explored in 
Baruteau & Lin (2010). Baruteau & Lin (2010) found that 



ip « 8.5 x 10~ 2 hs/cc 



(26) 



where the h term enters due to the mode lifetime being 
set by the speed of sound. Our model does not have a 
constant value for h, but we set this to 0.05 in this rela- 
tion. We use two disks in our simulations, with a values 
of 10~ 2 and 10" 3 , so following Baruteau & Lin (2010), we 
take ip 4.25 x 10" 4 and 1.34 x 10~ 4 , respectively. 

The sum of the n independent oscillations, <S> c ,m, spec- 
ifies the shape of the potential, <P = X^ c , m , that, in 
Laughlin et al. (1994), acted on the gas to generate den- 
sity perturbations, which in turn applied a force to the 
planet. We do not apply the potential to our gas disk, 
which is only evolved in one dimension. Instead, the 
potential is used to calculate the force acting on the 
planet directly. Ogihara et al. (2007) state that the force 
from such a turbulent disk can be expressed as 



turb 



-CV<2> 



(27) 



where C relates the size of the force acting from the po- 
tential onto the gas to the size of the force acting from 
the gas onto the planet 



64X> 2 



(28) 



The strength of the force scales with the surface den- 
sity of the gas at the planet's location, S(r). Eq. (27) can 
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FIG. 5. — Time histories of models H, which initially consist of planetary embryos with 0.75, 0.5, 1.0 and 1.25 M®, respectively. Each row depicts 
all four runs at a given time in its history-the top row represents the initial conditions, followed by 0.1, 0.2, 0.4, 1.5, and 2.5 Myr. Each panel depicts 
the semi-major axis and eccentricity of the bodies, with the area of a body scaling linearly with mass. Circles corresponding to masses of and 
lOMc, are shown for comparison in the bottom row. All but the 0.5 run H2 form a core of at least 10 Mp, by 1.5 Myr 



be used to calculate the force due to turbulent pertur- 
bations on a planet at position (r,8) at time t. Ogihara 
et al. (2007) showed that all modes with wavenumbers 
higher than 6 can be left out of the summation to deter- 
mine the potential <P, and so we follow this simplifica- 
tion, including only perturbations where 1 < m < 6 in 
our calculation of <P. 

2.5. N-BodyCode 

We simulate the behavior of multiple planetary em- 
bryos approaching a convergence zone, using the 
Bulirsch-Stoer N-body code described by Sandor et al. 
(2011), modified to include the additional forces from 
the gas disk on each body. The temperature and surface 
density profiles from the one-dimensional hydro simu- 
lations conducted previously are read into the N-body 
simulation at each timestep and used in the prescrip- 
tions described in Sections 2.2, 2.3, and 2.4 to determine 
the additional forces. 



The total force resulting from Type I migration, damp- 
ing, and turbulence are computed for each body at the 
begining of a Bulirsch-Stoer timestep, but are held con- 
stant in each refinement of the modified midpoint al- 
gorithm during a given Bulirsch-Stoer timestep. The 
Bulirsch-Stoer timesteps are a small fraction of the dy- 
namical time of the bodies, and are reduced during close 
encounters, so the simplification of holding these forces 
constant during a timestep does not change the overall 
behavior of the bodies. The net force acting on a body is 
given by 

F total = -Fnbody + ^typel + ^ damp + F W rb (29) 

where f n b dy represents the gravitational forces from 
the other bodies and the central star. The net force from 
the migration torques is 

Ftypel = t ^ PeI 0/ (30) 
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FIG. 6. — Time histories of runs with a Gaussian mass distribution centered at l.OMe (model Gl) , 0.75Me (G3), and 0.5M^ (G5) are presented, 
along with a superimposed distribution of the low density populations G2 & G4 (G8). The top row shows the initial conditions, followed by the 
system after 0.1, 0.2, 0.4, 1.5 and 2.5 Myrs. The results of core formation, and the time required to build a core, are comparable for the Gaussian 
and uniform mass distributions. 



while the components of F^ am p are given by Eqs. (22) 
and (23), and F tur b is derived in Eq. (27). 

A collision occurs if two bodies pass within a distance 
equal to the sum of the radii of the two bodies, where 
the radii are derived by assuming that all bodies have 
a bulk density of 2 gem -3 . Collisions result in the for- 
mation of one body with a mass and momentum equal 
to the sum of the two coliding bodies. This approxima- 
tion is appropriate for bodies with the masses of plan- 
etary embryos that have escape velocities substantially 
exceeding typical collision velocities. 

3. INITIAL CONDITIONS 

3.1. Gas Disk Parameters 

The evolution of the gas disk was described in 
Sect. 2.1. The parameters that set the initial profile of 
disk density and temperature are a, which sets the vis- 
cosity, the background temperature, Tj,, and the mass 



accretion rate, M. MHD simulations consistently sug- 
gest that a. has a value of 10~ 2 (e.g. Davis, Stone & Pes- 
sah 2010), though if the migration is taking place in a 
dead zone, the viscosity should be smaller, though non- 
zero (Fleming & Stone 2003; Oishi and Mac Low 2009) 
Our simulations use a disk with a = 10~ 2 , an accretion 
rate M = 1O" 7 M0 yr -1 , an external radius r ex t = 30 AU, 
and Tf, = 10K. The parameters for the gas disks used in 
our simulation are set so that the convergence zone lies 
within the simulation at all times. 

The photoionized wind has an escape radius r ? = 
5 AU for a 1 M Q star, and a photoevaporation rate set by 
a mass loss rate from the disk of M w = 3x 1O~ 8 M0 yr -1 . 
The accretion rate M, along with the photoevaporation 
rate M w and radius r„, determines the lifetime of the 
disk. Using the given values for our parameters yields 
a 0.08 M© disk with an 8 Myr lifetime. 
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FIG. 7. — Left panel. The formation history of model Sw3, which has 11 lMcg bodies, migrating through a population of 25 0.25M S bodies. The 
massive bodies do not merge with the lower-mass ones as the massive bodies rapidly migrate towards the convergence zone. Low mass bodies 
are briefly scattered, but then resume their slower migration rate. The massive bodies finish migrating to the convergence zone by 0.3 Myr, while 
the low-mass bodies continue to migrate until 1.0 Myr. An 8.5 Meg core is formed, though only two of the 25 0.25M^ bodies are consumed in 
forming it. Right panel. Extended history of the model. The swarm of smaller particles is seen to be long lived compared to the time required to 
form a planetary core. The swarm of low mass embryos does eventually merge with other bodies as they migrate inwards. 



3.2. Planet Embryos 

The location of the convergence zone, where F t y pe i — > 
0, is determined by an opacity transition that in turn is 
set by the background temperature and the viscosity of 
the disk. In a disk with Ty = 10 K and a = 10~ 2 , the con- 
vergence zone begins near 22 AU and moves inwards 
towards the star on the viscous evolution timescale. 
Since this is the region of interest, and the computation 
time using the Bulirsch-Stoer method scales as the num- 
ber of bodies squared, we distribute the initial semi- 
major orbital radii of planet embryos uniformly around 
the radius of the convergence zone. Most runs begin 
with planetary embryos in the region between 8 and 28 
AU. 

A consensus has not been reached on the initial mass 
function of planetary embryos. Starting with an initial 
population of 3000 equal mass (10 23 g) seed planetes- 
imals, and assuming a growth process dominated by 
binary collisions in the absence of gas drag and frag- 
mentation, Kokubo & Ida (1996) find a resulting mass 
distribution of power-law index 2.5 ± 0.4. The result 
of the opposite process, a collisional cascade from frag- 
mentation of a parent body or collection of parent bod- 
ies with subsequent self-similar grinding of the frag- 
ments, yields a quasi-steady mass distribution of power 
law 11/6 (Dohnanyi 1969; Tanaka et al. 1996). Lyra et 
al. (2008) found from hydrodynamical models of pro- 
toplanetary disks including interacting decimeter-sized 
particles and gas tides, that the resulting mass dis- 
tribution of lunar and Mars-sized embryos obeyed a 
power-law index of 2.3 ± 0.2, consistent with coagula- 
tion (Kokubo & Ida 1996) and, at 3c, with fragmentation 
(Dohnanyi 1969), suggesting a combination of both pro- 
cesses. However, none of these power laws match the 
observed mass distribution of the Asteroid and Kuiper 
belts (e.g., Cuzzi et al. 2010; Morbidelli et al. 2009), 
which are characterized by multiple power-laws (see fig 
1 of Morbidelli et al. 2009). Instead of adopting these 



power-laws, we use simple mass distributions, detailed 
below. We justify this by noting that the relaxed mass 
distribution is nearly independent of the initial mass 
distribution (Kokubo & Ida 1996). The parameters of 
the various models we ran are shown in Table 1. We 
detail them in this and in the next section. 

First, in our fiducial model (model HI in Table 1), the 
initial masses of the planetary embryos are all set to 0.75 
Mg . The innermost embryo is placed in a circular orbit 
with radius a\ = 8 AU, and each succesive embryo is 
placed so that neighbors i and j are initially separated 
by N x R m H, where the mutual Hill radius 

/ ntj + ntj \ 1/3 / fl; + ai \ 

r -«=(^mt) (V 1 )' (31) 

and cij = a-i + NR m H- In the fiducial model HI, we use 
N = 3. Several other runs were performed with varying 
masses of initial bodies, but maintaining the constraint 
of spacing consecutive bodies by three mutual Hill radii. 
These runs are all labeled H, for "Hill spacing", in Table 
1. 

Second, we used a Gaussian distribution of the initial 
masses of planet embryos. This is done by setting the 
mean mass and standard deviation for the distribution, 
with a floor of 0.1 M® and a ceiling of 2M®. These mod- 
els are labeled G, for "Gaussian", in Table 1. In these 
runs, the mutual Hill radii constraint is still used to de- 
termine initial locations, but the number of Hill radii 
between neighboring bodies is varied between runs to 
model low density and high density populations. Three 
further simulations are performed with an initial popu- 
lation of bodies resulting from superimposing two low 
density populations (models G7-G9). 

Third, we ran initial conditions consisting of two pop- 
ulations of planetary embryos with different masses, to 
investigate how massive, quickly migrating bodies in- 
teract with lower-mass bodies that migrate on much 
longer timescales. The low-mass population in this 
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FIG. 8. — Simulation histories of models Sw, where a population of planets evolves while interacting with a swarm of 25 low-mass (0.25 Mj) 
embryos. Column 1 shows the evolution of the swarm alone (model Swl). The more massive bodies quickly migrate to the convergence zone 
where they begin to form a planetary core. Only a small number of the 0.25Mc, bodies collide with anything, while most follow chaotic orbits in 
two distinct clusters, one interior and one exterior to the planetary core. These swarms eventually merge at later times, as seen in Fig. 7. 



model consists of 25 embryos with a mass of 0.25 
separated by 4 mutual Hill radii, while the massive pop- 
ulation consists of planets with masses of 0.75 or 1.0 
M9 separated by either 6 or 10 mutual Hill radii. The 
two populations are generated independently of one an- 
other, so a large and small body may initially be in close 
proximity to one another. These models are labeled Sw, 
for "swarm", in Table 1. 

Fourth, we use an initial configuration that scales with 
the amount of disk mass within the orbital radius of a 
planetary embryo. The mass of the planetary embryos 
is held constant, so the initial spacing between the semi- 
major axes of neighboring embryos is varied so that the 
mass of the gas in the annulus separating the two bod- 
ies is constant. Assuming that the dust-to-gas ratio in 
the disk is 0.01, and planetesimal formation turns 15% 
of the dust into planetary embryos, then the initial spac- 



ing between two neighboring bodies 

Ar= (6.67 x 10 2 )-^-. (32) 

where the numerical factor is the inverse of the product 
of the dust-to-gas ratio 0.01 and the planet formation ef- 
ficiency 0.15. These models are labeled SD, for "surface 
density", in Table 1. One simulation is performed us- 
ing a Gaussian mass distribution, spaced using the SD 
criteria. This simulation is included in the SD suite of 
simulations. 

The initial eccentricities and inclinations of planets 
are also selected randomly from Gaussian distributions. 
The inclination is the absolute value of a Gaussian dis- 
tribution with mean and 0.05° as the standard devi- 
ation.The mean value for the initial eccentricity is 0.05, 
with a standard deviation of 0.02. If the eccentricity is 
negative, a new value is generated until the value is 
non-negative. Due to the magnitude of the damping 
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FIG. 9. — The history of model SD1, a run with 0.75M^ planets ini- 
tially distributed such that the surface density of planets is 0.15% of 
the gas surface density. The results are comparable to the fiducial 
model, suggesting that mass distribution does not alter the qualita- 
tive description we have given of the formation of planetary cores. 

force, initial eccentricities and inclinations are found to 
be transient and do not play a significant role in the sim- 
ulations. 

4. RESULTS 

4.1. The Fiducial Model 

The history of the fiducial model (HI), with 23 bodies 
each having a mass of 0.75M(q and separated by 3R m n, 
is shown in Fig. 2. Initially, the motions of the embryos 
are dominated by the Type I torque, leading to rapid 
convergence. The convergence leads to close encoun- 
ters between neighboring bodies, resulting in collisions 
and orbit swapping. The close encounters and chaotic 
scattering and collisions persist until the number of re- 
maining planets is small enough to quiescently migrate 
inwards. 

The sharp peaks near 600,000 and 800,000 years are 
not planets being scattered out to large radii, but rather 
the result of solving for the orbital elements, assuming 
the planet is orbiting the central body, when in reality 
the bodies have formed a transient binary system. Such 
binary systems and satellite capture have been observed 
in other N-body simulations that include the effects of a 
gas disk, such as Cresswell & Nelson (2008). 

The fiducial model has a total initial mass of 17.25 
in planet embryos. In it, a 10. 5M® core forms at 
the zero-torque orbit in 1.29 Myr. No further collisions 
happen after this, as the remaining bodies are trapped 
in mean-motion resonances with the massive planetary 
core (Fig. 3). Two 3M® bodies are trapped in the 5:6 
resonance with the planetary core, and a 6M© body is 
trapped outside the orbit of the planetary core, in a 7:6 
resonance with it. An additional 3M® embryo orbits 
within the co-orbital region of the planetary core as a 
Trojan. Other simulations show similar behavior, with 
accretion eventually being halted by resonance trapping 
once the planetary core is massive enough. 

4.2. Variations of Planetary Embryo Mass 



The fiducial model begins with planetary embryos of 
0.75 M®. We also simulate disks populated with em- 
bryos of 0.5, 1.0, 1.25, and 1.5 Mg separated, as in the 
fiducial case, by 3R m u- These correspond to models H2- 
H5 in Table 1. Given the dependence of the Hill radius 
on the mass of the planet, systems consisting of more 
massive embryos contain a greater total mass. 

Since the magnitude of the migration torque, Eq. (12), 
varies as the square of the mass of the planet, and thus 
the acceleration varies linearly with mass, lighter em- 
bryos take longer to migrate to the convergence zone. 
This, combined with gravitational focusing, explains 
why the time required to build a core is a strong function 
of the initial mass of the embryos (Fig. 4). The fiducial 
model built a lOM® core in 1.28 Myr, while the run us- 
ing TOM® embryos reached a lOM© core in 0.43 Myr. 
The runs with 1.25 and 1.5 embryos formed cores 
in 0.21 and 0.15 Myr, respectively. The 0.5 M© embryos 
failed to produce a 10 M© body, but did form an 8 M® 
core after 2.22 Myr. 

To see if these results are characteristic of such a mass 
distribution, identical simulations to the ones above 
with initial planet masses of 0.75, 1.0, 1.25, and 1.5 M® 
are run, though the azimuthal distribution of the planets 
is varied, and so the chaotic interactions lead to different 
results. The dependence of the time required to assem- 
ble a massive core does hold for the five runs of each 
initial mass. There is, however, a wide upper bound for 
the time to form a core, so much so that each distribu- 
tion other than the 1.5 M© runs (model H5) had one run 
of the 5 fail to form a 10 M© core. These are not plotted 
in Fig. 4, but the resulting systems were interesting, of- 
ten consisting of two or more massive bodies coorbiting 
with one another stably for the rest of the 8 Myr simula- 
tion. 

The formation history of runs Hl-4 are shown in 
Fig. 5. The simulations with more massive embryos are 
seen to rapidly develop large planets, while the simula- 
tions consisting of embryos with initial masses less than 
1M® remain dispersed after 4 x 10 5 yr. The simulations 
demonstrate that, while the time required to form a core 
can vary by over an order of magnitude for embryos 
whose initial mass differs by a factor of two, the result- 
ing core and the arrangement of the remaining bodies is 
similar for all the simulations. The results presented in 
Fig. 5 are representative of the collection of simulations 
run with these initial conditions. 

Runs Gl-9 also suggest that the rate of core formation 
is highly dependent on the mass of the embryos. Two 
runs with a mean mass of l.OM© and a standard de- 
viation of 0.3M© (Gl and G2) built large cores within 
0.3 Myr. The initial masses for the two distributions 
are generated independently, which also varies the ra- 
dial distribution, since planets are initially separated by 
3R„,h- The core formation history of one of these runs 
is shown in Fig. 6, along with the formation histories of 
G3, G5, and G8. The results of these runs suggest that 
the trend seen in Fig. 4 can be used to estimate the time 
required to assemble a core from a Gaussian mass dis- 
tribution with a given mean mass. 

If the mean of the initial mass distribution is lowered 
to O.SM© (models G5 and G6) the time required to mi- 
grate bodies to the convergence zone increases signifi- 
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cantly, and when cores are built, they are smaller and 
take longer to form. The mass distribution is again 
Gaussian with a standard deviation of 0.2M®, though 
there are lower and upper mass limits of 0.1M© and 
2.OM0, respectively. Two runs were performed, with an 
independent mass distribution generated for each. The 
first run took 2.5 Myr to build a 6.9M© core, and is de- 
picted in Fig. 6 beside the history of a run composed of 
only 0.5M© planetary embryos. The second run took 3.2 
Myr to form a 9.3M© core. 

The runs with a population of 0.25M© bodies inter- 
spersed with more massive embryos did not result in 
the smaller bodies being swept up in collisions with the 
rapidly migrating large embryos. Instead, the larger 
embryos scatter the 0.25M© objects into high eccentric- 
ity orbits and continue to migrate towards the conver- 
gence zone where mergers between large bodies go on 
to create cores. This can be seen in Fig. 7 (left panel) by 
the two characteristic timescales for migration, which, 
aside from a few brief scattering events, occur more or 
less independently of one another. Following the for- 
mation of a core with mass ~ 5M®, smaller bodies are 
seen to collide with the core or become trapped as Tro- 
jans. Many, though, avoid close encounters with the 
planetary core and continue to interact chaotically with 
the other remaining small planetary embryos in two 
swarms, one interior and one exterior to the planetary 
core. Neglecting the (substantial) effect of gas accre- 
tion onto the planetary core, the model shows that these 
swarms eventually merge at later times, as seen in the 
right panel of Fig. 7, which depicts the whole 8 Myr of 
the simulation. Figure 8 depicts the evolution of the 
background swarm of 0.25M© by themselves (model 
Swl) as well as three cases where more massive bod- 
ies are interspersed into the low mass swarm (models 
Sw2-Sw4). The mass of the larger bodies, as well as the 
number of larger bodies, is not seen to have a notice- 
able effect on the migration and merging of the smaller 
bodies. 

4.3. Mass Distribution 

An alternative to determining the initial positions of 
planetary embryos by requiring them to be a set number 
of mutual Hill radii from their neighbor is to require that 
the distribution of mass in the form of planets mirror the 
mass distribution of the gas disk. This leads to a larger 
number of planets starting at radii exterior to the zero- 
torque orbit. In Fig. 9, we compare the history of fiducial 
run HI with run SD1, which differs from the fiducial run 
only in that the planetary mass is linearly proportional 
to the mass of gas within an annulus, following Eq. (32). 

The early period of convergence is skewed to larger 
radii in the latter, but as can be seen in Fig. 1, the torque 
exterior to the convergence zone is of the same magni- 
tude as the torque interior to the convergence zone. The 
time to migrate from initial locations to close enough 
to the convergence zone that N-body effects take over 
is not significantly affected by altering the initial radial 
distribution of planets. 

5. CAVEATS 

The prescription for Type I migration used in this pa- 
per was derived by Paardekooper et al. (2010) using 2D 



simulations. Those authors noted that the gas dynamics 
may act differently in 3D. 

Our simulations assume that the Type I migration 
torque derived for an isolated body in a non-isothermal 
disk also applies to many bodies that pass arbitrarily 
close to one another. This close proximity may signifi- 
cantly alter the migration torque or lead to a more rapid 
saturation of the torque. Such proximity is often short- 
lived, but in the case of a satellite being captured could 
lead to a change in the torque acting on the pair last- 
ing thousands of years. A preliminary two-dimensional 
hydrodynamical simulation of two bodies migrating in 
close proximity suggests that the wakes of the planets 
drive brief periods of large variations in the torque. This 
will be explored in future work. 

The initial convergence takes place while planets are 
separated from one another, so using the prescription 
for the torque on an isolated body is reasonable. Follow- 
ing convergence, close gravitational encounters domi- 
nate until a planetary core is formed. In several simula- 
tions, the resulting system had two large bodies that had 
a combined mass greater than 10M® but failed to pro- 
duce a single massive core. Studying such a system in a 
simulation that accounts for gap opening by the planets 
in the gas disk could shed light on how such an arrange- 
ment will evolve. 

Bitsch & Kley (2010) showed that the horseshoe 
torque is reduced as eccentricity is increased, and shuts 
down as the radial excursions associated with eccentric- 
ity exceed the width of the co-rotation region. Bitsch 
& Kley (2011) further showed that outward migration 
can be sustained for inclinations up to about 4.5°, after 
which migration stalls in general, owing to the lower 
densities away from the disk midplane. We assess the 
effect of eccentricity and inclination in the next section. 

In addition, the prescription for the co-rotational 
torque that we use considers only the unsaturated 
regime. We calculate for the parameters used in the 
present work that the effects of diffusion, explored in 
Paardekooper et al. (2011) for alpha disks, imply that 
only larger planets, with masses exceeding « 4M®, will 
experience sustained outward migration. Our choice of 
working with the unsaturated torque stems from the 
results of Baruteau & Lin (2010), Uribe et al. (2011), 
and Baruteau et al. (2011). These works measured the 
co-rotational torques in turbulent disks, finding that 
the stochastic turbulent fluctuations work toward keep- 
ing the co-rotational torque unsaturated even in locally 
isothermal simulations. Baruteau et al. (2011) have con- 
clusively shown that horseshoe dynamics exists in tur- 
bulent disks, and argue that it occurs when the ampli- 
tude of the U-turn drift rate exceeds that of the turbulent 
velocity fluctuations. This condition may or may not oc- 
cur for smaller planets, for which the U-turn drift rate is 
reduced. Unfortunately, neither Baruteau et al. (2011) 
nor Uribe et al. (2011) could numerically resolve the ex- 
pected horseshoe region width in this mass range. We 
thus recognize saturation for smaller planets as a possi- 
bility, but regard the question as unsettled. 

Our conclusions about resonance trapping neglect the 
effect of interactions with the wakes of the growing 
planetary core. These may be sufficiently strong to 
knock lower-mass bodies out of the resonances that we 
find. Even if the resonances do remain, our assumption 
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FIG. 10. — The eccentricity e normalized by the dimensionless half- 
width of the co-rotational region x s , for the protoplanets in four of our 
models. The eccentricity is pumped well above the critical value of 
e/x s = 1 associated with shutdown of the horseshoe drag (see text). 

about evolving the gas disk independently of the plan- 
ets breaks down once the planetary core has reached a 
mass near 10 M®. At that point it begins first to accrete 
gas, and then, when its Hill radius exceeds a disk scale 
height, to open a gap in the disk. Our models account 
for neither development, and so should be taken as only 
qualitatively indicative (e.g. the right panel in Figure 7) 
after any object reaches this mass. 

We note that the actual critical mass necessary for run- 
away gas accretion is a sensitive function of the lumi- 
nosity of the forming planet, which is determined by the 
release of gravitational binding energy of the accreted 
material — usually assumed to be provided steadily by 
planetesimals. In our model, this continuous luminos- 
ity is essentially zero, which would lead to a very small 
critical core mass (Papaloizou & Terquem 1999). The en- 
ergy is instead provided in discrete violent collisions. 
A model that incorporates a self -consistent treatment of 
gas accretion and gap formation by a growing planet 
with envelope luminosity provided by both planetesi- 
mal accretion and binary collisions would be required 
to fully assess the initial stages of giant planet forma- 
tion. 

5.1. Comparison to Recent Work 

The work presented in this paper examines the same 
regime as that studied by HN12, which was made pub- 
lic after the initial submission of our work. Although 
both works deal with the impact of co-rotation torques 
and convergent migration on the growth of protoplan- 
ets, there are differences in the models that are worth 
highlighting, in order to understand the different results 
obtained. 

First, we discuss the gas evolution model. In the 
present work we model the gas density by equating the 



Fig. 11. — The inclinations i normalized by the aspect ratio h of the 
disk, for the same models as in Fig. 10. The inclinations are pumped 
well above the critical value of i a n = 0.08 (i.e., i CT n/h = 1.6) beyond 
which migration stalls (see text). This is a potentially important phe- 
nomenon, perhaps countering the shutdown of co-rotational torques 
due to finite eccentricity shown in Fig. 10. 



radial velocity with the viscous inflow, the azimuthal 
velocity with the Keplerian value, and solving the re- 
sulting diffusion equation in the presence of photoevap- 
oration. The initial condition for density is set by defin- 
ing a global value for mass accretion rate and viscosity 
(a) value, which in turn define the surface density (Pa- 
paloizou & Terquem 1999). Having defined the density, 
the disk is assumed to maintain thermal equilibrium be- 
tween heating by viscous dissipation and cooling via 
black body radiation in a background radiative field of 
temperature Tj,, according to equation (5). In HN12, 
the density and temperatures are set as global power- 
laws (that may drop well below the physical value of 
Tf, ~ 10 K set by the surrounding molecular cloud, and 
even below the temperature of the cosmic microwave 
background). The viscous-photoevaporative evolution 
of the density is mimicked by an exponential decay with 
a fixed e-folding time Tdisk/ me temperature is evolved 
dynamically assuming thermal diffusion in the radial 
direction. 

Next we consider the growth of the protoplanets. 
While our work concerns the formation of a core of a gi- 
ant planet in the vicinity of a convergent migration zone 
via binary collisions, HN12 attempt a more comprehen- 
sive approach. In addition to collisions with planets 
of terrestrial size, their protoplanets also accrete from a 
disk of 10 km planetesimals subject to Stokes drag. They 
also include an approximate scheme for gas accretion 
based on fits to the fast ID giant planet formation model 
of Movshovitz et al. (2010). In their model, gas cap- 
ture is a function only of core mass and time, starting 
at roughly 3 M®, with runaway gas accretion setting in 
at 30 M®. The gas accretion is artificially cut off at either 
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the Jovian mass or at the gas isolation mass, whichever 
is lower. 

A minor difference between our work and HN12 lies 
in the inclination and eccentricity damping. Both works 
use prescriptions that are based on the analytical formu- 
lae derived by Tanaka et al. (2002) and Tanaka & Ward 
(2004). HN12 use the same inclination damping, but 
a smoother version for eccentricity damping. We use 
the model of Cresswell & Nelson (2008) that improves 
on Tanaka & Ward (2004) by including empirical fits to 
the transition from exponential to quadratic damping at 
higher eccentricities and inclinations. 

Finally we consider the prescription for migration of 
the protoplanets. HN12 model saturation effects, while 
we chose to work with the fully unsaturated torque, as 
explained in the previous subsection. Another major 
difference is our inclusion of stochastic migration (Nel- 
son & Papaloizou 2004; Nelson 2005) caused by the in- 
clusion of turbulent jitter (Sect. 2.4). 

As a result, HN12 find a bimodal formation process. 
Cores that grow too massive too fast migrate into the 
star, whereas cores of moderate mass slowly migrate 
outward to large distances, accreting gas and becom- 
ing giant planets at wide orbits, mostly beyond 10 AU. 
Super-Earths or Neptune-mass planets were not easily 
formed. 

The difference in the results that we find are readily 
understood in terms of these differences in the mod- 
eling. First, without a background radiation field of 
temperature T b below which the disk cannot cool, the 
model of HN12 lacks the transition to isothermality that 
gives rise to a limit of outward migration in our models, 
usually between 10-20 AU, depending on the initial pa- 
rameters used. Second, their inclusion of saturation in 
HN12 introduces a mass-dependency to both the loca- 
tion and existence of convergence zones. Rapid inward 
migration occurs for planets outside of a narrow range 
in mass. Lastly, our inclusion of turbulent jitter dynami- 
cally excites the resonant convoy, leading to more grav- 
itational interaction than seen in HN12, and faster core 
growth. These differences highlight the fact that planet 
migration and growth depend very sensitively on the 
disk model adopted. 

In two out of their suite of forty simulations, HN12 
apply a eccentricity cutoff to the co-rotational torque, 
taking in to account the findings of Bitsch & Kley 
(2010). They find that the eccentricity excursions of the 
protoplanets easily exceed the threshold of e/x s = 1, 

where x s = 1.1/ j l ^^/qjh is the dimensionless half- 
width of the co-rotational zone, therefore quenching the 
co-rotational torques. We measure e/x s in four of our 
models. The result is shown in Fig. 10. We see that 
the eccentricies excited by the N-body interactions and 
the turbulent jitter on the protoplanets far exceeds the 
stringent threshold of e/x s — 1, confirming the result 
of HN12. However, the same processes that account 
for eccentricity pumping should also lead to inclina- 
tion growth. As shown by Bitsch & Kley (2011), migra- 
tion stalls in general when the inclination rises above 
a threshold z crit « 4.5° « 0.08. We plot in Fig. 11 i/h 
for the same models as in Fig. 10. A comparison of the 
two figures suggests that eccentricity and inclination are 
indeed correlated. The aspect ratio being h = 0.05, we 
have JcritA = 1.6. Figure 11 shows that the inclinations 



grow well beyond this modest threshold. This indicates 
that although the high eccentricities would lead to a 
shutdown of the co-rotational torques and resuming of 
fast inward migration, the associated high inclinations 
would mitigate the effect by also shutting down the 
Lindblad torque and, consequently, migration in gen- 
eral. 

6. CONCLUSIONS 

We have studied the evolution of multiple planets 
migrating in a non-isothermal gas disk, assuming un- 
saturated migration torques (Paardekooper et al. 2010). 
Initially, the planets migrate towards the convergence 
zone as was described by LPM10 for individual planets 
in such a disk. Once the planets are near the conver- 
gence zone, however, the close proximity to other plan- 
ets leads to chaotic interactions and collisions, even- 
tually resulting in the formation of several large bod- 
ies that often merge to form a single massive plane- 
tary core. The planetary core dominates the dynamics 
of the remnant bodies, and the smaller bodies tend to 
get trapped in mean-motion resonances with the core. 
Some become trapped as Trojans with the core, and the 
whole system migrates inward on the viscous timescale 
as the disk evolves, maintaining the resonances. As the 
goal of this investigation was to study planetary em- 
bryos undergoing Type I migration, the simulation does 
not model accretion onto the planetary core or Type II 
migration. Both should be treated in future work. 

Since the migration torque scales with the mass of the 
planet, simulations that initially consist of larger mass 
planetary embryos produce planetary cores rapidly, 
while simulations that start with less massive bodies of- 
ten take more than 1 Myr to form a core. The rate of core 
formation does not appear to depend on the initial dis- 
tribution of planets in the disk, but does depend on the 
mass of the planetary embryos that coalesce to form the 
core. A realistic protoplanetary disk would consist of 
planetary embryos with a range of masses, but we have 
shown that even in the case of a Gaussian distribution 
of embryo masses the trend of large embryos migrating 
faster and merging earlier holds. 

The location of the convergence zone is determined by 
the temperature profile of the gas disk. This, in turn, is 
determined by the turbulent viscosity as parameterized 
by a, the total disk mass, and background illumination 
from the central star and neighboring stars. The conver- 
gence zone falls at larger radii for larger values of a and 
the disk temperature. More massive stars have hotter 
disks, leading to convergence zones at larger radii. Di- 
rect imaging has indeed revealed a population of giant 
planets at radii of order 100 AU (e.g. Marois et al. 2008; 
Kalas et al. 2008; Oppenheimer et al. 2008; Marois et al. 
2010) around stars significantly more massive than the 
Sun. Our results offer a possible means of building gas 
giant planets via fast core accretion at such large radii. 
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